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ABSTRACT 

We use ZEUS-MP to perform high resolution, three-dimensional, super- 
Alfvenic turbulent simulations in order to investigate the role of magnetic fields 
in self-gravitating core formation within turbulent molecular clouds. Statistical 
properties of our super-Alfvenic model without gravity agree with previous sim- 
ilar studies. Including self-gravity, our models give the following results. They 
are consistent with the turbulent fragmentation prediction of the core mass dis- 
tribution of Padoan & Nordlund. They also confirm that local gravitational 
collapse is not prevented by magnetohydrodynamic waves driven by turbulent 
flows, even when the turbulent Jeans mass exceeds the mass in the simulation 
volume. Comparison of results between 256 3 and 512 3 zone simulations reveals 
convergence in the collapse rate. Analysis of self-gravitating cores formed in the 
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simulation shows that: (1) All cores formed are magnetically supercritical by at 
least an order of magnitude. (2) A power law relation between central magnetic 

1 /2 

field strength and density B c oc p c is observed despite the cores being strongly 
supercritical. (3) Specific angular momentum j oc i? 3//2 for cores with radius R. 
(4) most cores are prolate and triaxial in shape, in agreement with the results 
of Gammie et al. We find a weak correlation between the minor axis of the core 
and the local magnetic field in our simulation at late times, different from the 
uncorrelated results reported by Gammie et al. The core shape analysis and 
the power law relationship between core mass and radius M oc i? 2 - 75 suggest the 
formation of some highly flattened cores. We identified twelve cloud cores with 
disk-like appearance at a later stage of our high-resolution simulation. Instead of 
being tidally truncated or disrupted, the core disks survive and flourish while un- 
dergoing strong interactions. We discuss the physical properties of these disk-like 
cores under the constraints of resolution limits. 

Subject headings: ISM: clouds — ISM: kinematics and dynamics — ISM: magnetic 
fields — stars: formation — turbulence — MHD — methods: numerical 



1. Introduction 

In the standard theory of isolated star formation, the presence of magnetic fields in 
molecular clouds plays a vital role. Based on the Jeans argument, observed high-density 
molecular clouds should collapse within a few freefall times to form stars. Instead, the 
molecular clouds were thought to be quasi-stable. Turbulence in molecular clouds and mag- 
netic field support were proposed to explain these apparently stable clouds (e.g. McKee 1999; 
Williams, Blitz, & McKee 2000; Vazquez-Semadeni et al. 2000). Molecular clouds were sug- 
gested to be supported magnetostatically (e.g. Mouschovias & Spitzer 1976; Fiege & Pudritz 
1999) or dynamically by magnetohydrodynamic (MHD) waves, especially Alfven waves (e.g. 
Dewar 1970; Shu, Adams, & Lizano 1987). Alfven waves are transverse waves and produce 
perturbation perpendicular to the mean magnetic field. McKee & Zweibel (1995) pointed out 
that Alfven waves can lead to an isotropic turbulence pressure to counteract gravitational 
collapse, provided that the waves are neither damped nor driven. However, this mechanism 
requires a negative radial gradient in wave sources in the cloud in order to support the cloud 
from collapsing (Shu, Adams, & Lizano 1987). 

Observationally, the question of whether magnetic fields are sufficiently strong to sup- 
port molecular clouds alone from collapsing remains unresolved. Crutcher (1999) concludes 
that static magnetic fields are insufficient to support the observed clouds and that MHD 
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waves are equally important in cloud energetics. McKee (1999) however points out that 
ambipolar diffusion might already have altered the mass-to-ffux ratio observed by Crutcher. 
Observations by Bourke et al. (2001) yield similar conclusions to Crutcher (1999), although 
they caution that the measured magnetic field strength may be significantly lower than the 
true values due to the low volume filling fraction of dense cores. Nakano (1998) argues 
that observed cloud cores cannot be magnetically supported. A magnetically supported core 
should not have a column density much higher than its surroundings, and the core can- 
not maintain a large velocity dispersion inside it for its whole lifetime. However, observed 
cloud cores generally have superthermal velocity dispersions, and column densities much 
higher than their surroundings, suggesting that most cloud cores are not magnetostatically 
supported. 

Significant progress in the understanding of the roles of turbulence and MHD waves in 
supporting gravitationally unstable cloud cores has been made recently using 3D numerical 
simulations (e.g. Vazquez-Semadeni, Passot, & Pouquet 1996; Mac Low et al. 1998; Padoan 
& Nordlund 1999a; Ostriker, Gammie, & Stone 1999, for a review, see Mac Low & Klessen 
2004). Klessen, Heitsch, & Mac Low (2000, hereafter Paper I), using models computed with 
both ZEUS-3D and a smoothed particle hydrodynamics (SPH) code, conclude that hydrody- 
namical turbulence can prevent global collapse but not local collapse under typical molecular 
cloud conditions. Heitsch, Mac Low, & Klessen (2001, hereafter Paper II) include magnetic 
fields and conclude that local collapse cannot be prevented much longer than a global free-fall 
time by magnetized turbulence, in the absence of mean-field support. However, they find that 
magnetic fields delay local collapse by decreasing local density enhancements via magnetic 
pressure. Stars begin to form quickly when local density enhancements collapse. This result 
favors the dynamical picture of molecular clouds being a transient feature in the interstellar 
medium (Ballesteros-Paredes et al. 1999; Elmegreen 2000; Hartmann, Ballesteros-Paredes & 
Bergin 2001) instead of existing for many free-fall times. 

Another interesting result from the 256 3 resolution simulations of Papers I and II is that 
several cores begin to evolve into flattened (or disk-like) objects, which leads to the specula- 
tion that they were beginning to resolve the formation of protostellar cores. Unfortunately, 
the resolution and number of disk-like cores were not high enough to perform a statistically 
meaningful study. Therefore, these flattened cores were not discussed in Paper II. 

Simulations suggest that supersonic turbulent flows, driven at large scale by supernovae 
and density waves, produce shock-compressed gas sheets or filaments that fragment into 
dense cores, and drive the observed supersonic turbulence in the clouds (Mac Low & Klessen 
2004). Differential rotation of galactic disks (Fleck 1981) provides another plausible driving 
mechanism to maintain interstellar turbulence, especially at galactic radii with little mas- 
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sive star formation. The coupling from large-scale shear down to turbulent scales could be 
through the magnetorotational instability (e.g. Balbus & Hawley 1998; Sellwood & Balbus 
1999; Kim, Ostriker, & Stone 2003). Ionizing radiation (McKee 1989; Vazquez-Semadeni, 
Passot, k Pouquet 1995; Bertoldi & McKee 1997; Kritsuk & Norman 2002), and stellar 
winds from massive stars may also be important, but supernovae (e.g. Norman & Ferrara 
1996; Avillez 2000) probably provide most of the energy required to maintain interstellar 
turbulence as well as the turbulence within molecular clouds (Mac Low & Klessen 2004). 

Recent numerical studies indicate that supersonic turbulence may play an important role 
in shaping the initial mass function (IMF) of the stellar population that forms in molecular 
clouds (Padoan & Nordlund 2002). Observations of the IMF appear consistent with a general 
form having a power law in the high-mass wing and then flattening and turning over at the 
low mass end (e.g. Kroupa 2002, and references therein). This is an important constraint on 
star formation theories. Interestingly, and perhaps significantly, the probability distribution 
function (PDF) of gas density in isothermal supersonic turbulence is well approximated by a 
log normal distribution. Padoan, Nordlund, & Jones (1997) and Padoan & Nordlund (2002) 
suggest that the IMF and density PDF in turbulent clouds are intimately related, and predict 
that the mass distribution of dense cores is controlled by the density, the Mach number, and 
velocity dispersion of the supersonic turbulence. Previous simulations of molecular cloud 
collapse (e.g. Padoan et al. 2001; Klessen 2001; Bate, Bonnell, & Bromm 2003; Gammie 
et al. 2003) indicate that the mass distribution of cores (or clumps) roughly resembles a 
log normal distribution, but with the high-mass wing following a power law approximately 
consistent with the Salpeter law. However, these simulations either lack magnetic fields or do 
not have a sufficient number of gravitationally bound cores for accurate statistical analysis, 
because of low numerical resolution. Higher resolution MHD turbulence simulations are 
required to verify the importance of super-Alfvenic turbulent fragmentation to the stellar 
IMF. 

In this paper, we report on results of our latest three-dimensional (3D), turbulent sim- 
ulation of a magnetized molecular cloud, increasing the spatial resolution in each dimension 
by a factor of two relative to the simulations in Paper II, using the ZEUS-MP code. The ex- 
cellent parallel scalability of ZEUS-MP allows us to study MHD turbulent molecular clouds 
on a 512 3 grid, sufficient to study the global physical properties of cores and their formation 
process. In §2, we describe the techniques, assumptions, and parameters used in our sim- 
ulations. In §3, we discuss the statistical properties of the super-Alfvenic turbulent cloud 
before gravitation is turned on. The global physical properties of self-gravitating cores and 
their time evolution are reported in §4. We also examine the turbulent fragmentation theory 
by comparing its predictions to our simulation results. In §5, we perform a resolution study 
and investigate the physical properties of some disk-like cores observed at later stages of our 
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simulation. We discuss our results and present conclusions in §6. 

2. Simulations 

Our simulations are performed on 256 SGI/Cray Origin 2000 processors at the National 
Center for Supercomputing Applications using the ZEUS-MP code. ZEUS-MP is a multi- 
physics, massively-parallel, message-passing code for astrophysical fluid dynamic simulations 
in 3D, developed by the Laboratory for Computational Astrophysics (lca.ncsa.uiuc.edu) at 
NCSA. ZEUS-MP is a distributed memory version of the shared-memory code ZEUS-3D, 
which uses block domain decomposition to achieve scalable parallelism (Norman 2000). The 
code includes ideal hydrodynamics, ideal MHD, and self-gravity. Self-gravity is implemented 
using the FFTW library (Frigo & Johnson 1998). 

Our scale-free MHD turbulent simulations are basically the same as model Ehlw in 
Paper II, but performed at resolutions up to 512 3 . All parameters are given in normalized 
units, where physical constants are scaled to unity. The total mass in the box M = 1 and the 
side length L = 2. The average density p Q = 0.125. The initial conditions of our simulation 
are identical to model Ehlw in Paper II. In this model, the sound speed c s = 0.1. The 
corresponding Jeans length Aj = 0.5, and Jeans mass 

M J = PoAj = (-^y^ Po 1/2 °l = 0-0156. (1) 

Therefore, the mass in the box M = 64 Mj and the length of the box L = 4 Aj. The 
magnetic field is initially uniform and along the z-direction with strength B Q = 0.188. The 
ratio M/M CT = 8.3, which means our entire cloud is magnetically supercritical (see e.g. 
Mouschovias & Spitzer 1976). The ratio of thermal to magnetic pressure is 0.9 and the rms 
Mach number for this parameter set is M rms = 10. Isothermality is assumed throughout the 
simulation, which is justified as the cooling time is much shorter than the dynamical time 
inside high density molecular clouds. No ambipolar diffusion is included in the calculation. 

Periodic boundary conditions are applied in all three dimensions. In effect we simulate 
a small portion of a larger cloud. We do this for simplicity, to avoid addressing the structure 
of the non-isothermal boundaries of molecular clouds. We do not expect that our results will 
be contradicted by unified models that include larger-scale flows in the interstellar medium. 
In our analysis, we examine column-densities through one realization of our periodic cube, 
which should give statistically valid answers to the questions we ask for molecular clouds of 
finite size. 

The procedure of setting up the turbulent flow is explained in Mac Low (1999). We 
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briefly review it here. We assume that each component of the velocity perturbation is a 
Gaussian random field with flat power spectrum in the range 1 < < 2. Random phases 
and amplitudes are generated in that spherical shell in Fourier space, and then transformed 
back into real space to generate each component of the driving velocity perturbation. When 
all three components are obtained, the amplitude of the velocity is scaled to a desired initial 
root mean square velocity, which is unity in our case. This velocity field is then normalized at 
each time step to ensure a constant kinetic energy input rate and added in as a perturbation 
to the evolved velocity field. The dynamical behavior of isothermal self-gravitating gas is 
scale free. Equations (2) to (4) can be used to scale back to astrophysical units with a mass 
scale of the thermal Jeans mass, length scale given by the thermal Jeans length, and a time 
scale given by the free-fall timescale (see Paper I), 

L = 0.89pc ( ,° a - ) ( A n - ) ~ 1/2 (2) 
P VO^kms- 1 / V10 4 cm- 3 / V ; 

M = 413M ( ^-—] 3 ( A U - ) ^ (3) 

° VO^kms-V Vl0 4 cm- 3 / V ; 

Therefore, if we choose a high-density scaling such as the BN region in Orion, in which 
c s = 0.2 km s" 1 and the number density n = 10 5 cm -3 , the region in the simulation box is 
0.28 pc in size. The total mass in this region is 130.6 M and one system time unit, which 
is 0.65 free-fall time, is 0.07 Myr. 

A total of 110,000 CPU hours were required for the 512 3 zone simulation to run to 4.325 
time units, or ~ 2.82 r ff . In our case, r ff = 1.53 system time units. At the beginning of the 
simulation, no gravitational forces are applied. Turbulence is allowed to develop from the 
relatively smooth initial conditions until t = 2.0 time units. Gravitation is then turned on. 
Density fluctuations generated by the supersonic turbulence in converging and interacting 
shock fronts that locally exceed the Jeans limit begin to contract. In the 256 3 simulations 
reported in Paper II, there were only a dozen gravitationally collapsing cores, which were 
under-resolved and insufficient for statistical analysis. In our 512 3 simulation, there are 15 
cores that satisfy our definition of gravitationally bound at t — 2.0 when gravity is turned 
on. By the end of the simulation at t — 4.325, we have 83 gravitationally bound cores. The 
simulation was terminated then because too many cells violate the Jeans condition (Truelove 
et al. 1997) which dictates that the local Jeans length must be resolved by several zones. 
Also, the densities of some cores became so high that the isothermal assumption would be 
invalid. Henceforth, we adopt the same time convention as in Paper II, where the time t — 
is defined as the time that gravity is turned on. According to this definition, the simulation 
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is terminated at t — 2.325 system time, or 1.52 r ff . We will use free-fall time as the time 
unit for the rest of the paper. 

3. Statistical Properties 

3.1. Power Spectra 

Kolmogorov (1941) heuristically derived a simple scaling for the ID energy spectrum 
in incompressible turbulence, E{k) oc A;" 5 / 3 , where k is the magnitude of the wavenumber. 
Astrophysical fluids are magnetized and highly compressible. Dynamically important mag- 
netic fields and strong shocks should interfere with eddy motions, which in turn will affect 
the final velocity power spectrum. However, observations and some simulations (e.g. Pa- 
per I; Lazarian, Pogosyan, & Esquivel 2002; Cho, Lazarian, & Vishniac 2003) suggest that 
MHD compressible turbulence also closely follows the Kolmogorov spectrum. The theory of 
Goldreich & Sridhar (1995) points out that it is the perpendicular motion of sub-Alfvenic 
incompressible turbulence that behaves as Kolmogorov, E(k) oc £r 5/3 . Some other MHD 
turbulence simulations suggest that the scaling law for MHD compressible turbulence is 
not exactly Kolmogorov. Boldyrev, Nordlund, & Padoan (2002) reported that the power 
spectrum of the solenoidal component of the velocity field in a super-Alfvenic compressible 
turbulence is E{k) ~ A;~ L74 , steeper than the Kolmogorov spectrum. In the latest simula- 
tions by Vestuto, Ostriker, & Stone (2003) on 256 3 and 512 3 grids, the slope of the power 
spectra varies from Kolmogorov to even larger than the E(k) ~ k~ 2 Burgers' model of 
shock-dominated turbulence (Burgers 1974) for models ranging from sub-Alfvenic to super- 
Alfvenic. The measurement of the slopes in their simulations might be affected by the short 
(or nonexistent!) inertial range of their power spectra, however. 

Figure la shows the power spectra of density, velocity, and the energies of our turbulent 
cube at t — 0, after the turbulence has fully developed, but before gravity has been turned on. 
The long dashed line is the driving spectrum. The power spectrum of velocity averaged over 
all directions is basically a Kolmogorov power spectrum (the thin straight line in Figure la). 
The inertial range of the velocity power spectrum spans more than an order of magnitude 
of the wavenumber (2 < \k\ < 30). The power spectrum of the magnetic field shows a 
similar result. The density spectrum has a shallower slope with index ~ -1.15. Both the 
kinetic and potential energies have shallower slopes because of the modulation by the density 
spectrum. Figure lb is a compensated plot of the velocity power spectrum multiplied by 
A; 5 / 3 that highlights the excellent fit of the Kolmogorov power spectrum to the velocity power 
spectrum. We also plot the compensated velocity power spectrum for the 256 3 resolution 
simulation for comparison. The inertial range of the turbulent flow in this lower resolution 
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simulation is shorter but still close to the Kolmogorov power spectrum. 

3.2. Density PDF and Core Mass Spectrum 

The PDF provides important and complementary information to the power spectrum. 
The PDF of isothermal turbulent flows can be approximated by a log normal function (e.g. 
Vazquez-Semadeni 1994; Padoan, Nordlund, & Jones 1997; Scalo et al. 1998; Nordlund & 
Padoan 1999). With a log normal PDF of mass density, most of the material concentrates 
in a small fraction of the total volume in the simulation. Figure 2 shows the best log 
normal fit to the PDFs of all the cells in the 256 3 and 512 3 simulations at t — 0. The 
standard deviation of the best fit is 1.46 ± 0.06 and 1.42 ± 0.03, as defined in equation (6) 
in Padoan & Nordlund (2002). Using equation (9) in Padoan & Nordlund (2002), a 2 = ln(l 
+ Mj^ rms rj 2 ), we calculate r\ to be 0.41 and 0.38, respectively, a little smaller than the value 
i] ~ 0.5 reported in Padoan & Nordlund (2002) from numerical experiments. 

Figure 3 shows the projected density (column density) images of the simulation, in 
logarithmic scale, along the x-, y-, and z-direction, which shows the dumpiness expected 
from a log normal PDF. The first row of plots are at t — 0, just before gravitation is 
switched on. The second and third rows are at t = 0.65 r ff and t = 1.31 r ff , respectively. 
At t — (Figures 3a-c), filamentary structures dominate the view and we see some higher 
density nodes on the filaments. Later (Figures 3d-f), gravity breaks the filaments into pieces 
and higher density clumps slowly increase their mass by accretion. The peak density at this 
moment is about 1000 times the average density. At t = 1.31r ff (Figures 3g-i), many cores 
with high column-densities have formed. The peak density is about four orders of magnitude 
above the average density and the total mass in cores is less than 9% of the total mass in 
the box. Qualitatively, this picture does not change during the subsequent simulation time 
(up to t = 1.52r ff ), except that more cores form and their peak densities grow. Most of the 
volume in the simulation corresponds to low density voids. The dumpiness can also be seen 
in Figure 4 which shows a volume rendering of the logarithmic density and the magnetic 
pressure of the turbulent box in 3D at t — 1.31 tr. An animation of the evolution of density 
and magnetic pressure of the 512 3 simulation is included in the electronic version of this 
paper. 

Padoan & Nordlund (2002) propose that the stellar IMF is a direct result of turbulent 
fragmentation of molecular clouds. As shown in Figure 3, the gas is compressed into ID 
filaments and 2D sheets by supersonic turbulence. Therefore, the core mass distribution 
will be dependent on the jump conditions for isothermal shocks in a magnetized gas. The 
typical size of a dense core is comparable to the thickness of the post shock gas. Based on 
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this, Padoan & Nordlund (2002) calculate by taking into account the scale-dependent Mach 
number that the mass distribution of cores formed inside a supersonic MHD turbulent cloud 
is 

N(m)dlog(m) oc m~ 3 ^ 4 ~^ dlog{m) (5) 

where (3 is the spectral power index. With our 512 3 simulation, we have sufficient resolution 
to check the turbulent fragmentation prediction. We proceed as follows. 

The clumpy structure inside the turbulent molecular cloud in our simulation is deter- 
mined using the algorithm CLUMPFIND described by Williams, De Geus, & Blitz (1994). 
The algorithm is also described in Paper I. A "core" is defined as a gravitationally bound 
region. We define it to consist of a set of connected zones with average density exceeding 
the mean density expected for isothermal shocks, p > M r 2 ms p , potential energy exceeding 
kinetic energy, \Epf t l \ > E^, and core mass exceeding the local Jeans mass, M core > Mj(p). 
We use logarithmic density contours in CLUMPFIND to get a wide enough density range, 
as in Paper I. We believe the core definition chosen here is more rigorous than the simple 
gravitationally unstable overdensity criterion of Padoan et al. (2001). However, our defini- 
tion identifies fewer cores. This definition of core will be applied to all results on cores in this 
paper. We also define "clumps" to be overdense regions that do not satisfy the other criteria. 
We choose the overdensity threshold for clumps somewhat arbitrarily to be p > 0.1M r 2 ms p o 

In Figure 5, we plot the mass spectra of cores (dot-dashed line) and clumps (dashed 
line) in the 256 3 simulation and the spectrum of cores (solid line) in 512 3 simulation at t — 
0, normalized by the Jeans mass Mj(p D ). There are only 11 cores in the 256 3 simulation. 
Therefore, we also plot the clump spectrum for comparison. Many of these clumps could be 
just density fluctuation created by the isothermal shocks that will be destroyed later. The 
256 3 clump spectrum and 512 3 core spectrum both show a well-defined "universal" IMF (see 
the latest review by Kroupa 2002). From Figures la and lb, the velocity power spectral 
index (5 ~ -5/3. Substituting (3 into equation (5), we expect the slope of the high-mass wing 
of the core mass spectrum will be -1.29. This spectral index is plotted as thin line in Figure 5 
and is consistent with the high-density wings of all the spectra, as predicted by the turbulent 
fragmentation. Even the 256 3 core spectrum appears consistent with the prediction. Note 
that the number of cores in Figure 5 is still too few to have an accurate fitting unless we 
relax the definitions of the cores to clumps. 
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4. Cloud Core Evolution 

4.1. Delay of Core Collapse 

Paper I and Paper II concluded that local collapse could not be prevented even when 
the turbulent velocity field carries enough energy to balance the gravitational contraction 
on global scales. Paper II found that the comparison between 128 3 and 256 3 MHD results 
did not show full convergence, with increasing resolution delaying local collapse (see Figure 
6 in Paper II). This occurs because increasing the resolution produces thinner shocks with 
peak densities closer to the theoretical value. The resulting deeper potential wells accrete 
more mass and collapse more quickly. At the same time, the simulation can better follow the 
short-wavelength MHD Alfven waves, which could create an isotropic pressure to counteract 
the gravitational collapse. From the 256 3 simulations in Paper II, it appears that the MHD 
waves are able to delay the collapse but not to prevent it. In our new 512 3 simulation, 
this imbalance between gravitation and Alfven wave pressure is verified. Figure 6 combines 
our new results with the previous 64 3 , 128 3 , and 256 3 results. We plot the collapsed mass 
M*, defined as the sum of the core masses, versus time. Collapse still occurs in all cases. 
Comparison between the 256 3 and 512 3 models suggests that the 256 3 result of Paper II is, 
in fact, converged. Remaining discrepancies are produced by the chaotic nature of the flow, 
as was shown by comparisons of different realizations of the same resolution model in Paper 
I. 



4.2. Core Formation and Evolution 

In the models, we find that self-gravitating cores form out of initially overdense clumps 
created by supersonic turbulence. We also find that dense clumps can be destroyed by 
passing shocks as well as by the merging of clumps to form larger and denser objects. The 
formation and destruction of cores is thus a complex dynamical process. The core formation 
and destruction process is shown in Figure 7 in terms of number of cores formed as a function 
of time. Instead of monotonically increasing, the core number drops sometimes because of 
merging or destruction. 

Figures 8 to 10 show the evolution of the physical properties of cores after gravitation 
is turned on and before the Jeans criterion Truelove et al. (1997) is violated. 

Figure 8a and 8b show the ratio of core mass to the mass that the local magnetic field 
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can support, M/M CT . The ratio 
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where iV is the column density; the geometrical correction factor F = 2 for a uniform sphere 
and 3 for an isothermal disk. B\ os is the line of sight magnetic field strength; G is the 
gravitational constant; and c$ ~ 0.12 for a spherical cloud or 0.16 for an isothermal disk 
(Bourke et al. 2001). In these two plots, we choose F = 3 and c$ = 0.16 for the extreme case. 
The ratios exceed an order of magnitude at all time. If we assumed spherical cloud geometry, 
the ratios would be even higher. Therefore, all the cores are magnetically supercritical during 
the gravitationally collapsing phase, despite the action of gravitational fragmentation, which 
could be expected to reduce the mass to flux ratio (Mestel & Spitzer 1956). 

Figures 8c and 8d show the relationship between the central magnetic field strength 
B c and the central core density p c at t = and t = 0.65 r ff respectively. The power law 
of index 0.5 is also plotted. Interpretations of B c oc pj 2 can be that: 1) cores are mag- 
netically subcritical (e.g. Mouschovias 1991; Basu & Mouschovias 1994, 1995), 2) cores are 
magnetically supercritical but the Alfvenic speed is about constant (~ 1) for the cores (e.g. 
Bertoldi & McKee 1992; Crutcher 1999), 3) cores are magnetically supercritical and dense 
cores tend to accrete mass along magnetic field lines and reduce their magnetic flux to mass 
ratio efficiently, even in the absence of ambipolar diffusion (Padoan & Nordlund 1999b). 
The cores in the simulations are all magnetically supercritical and there is no ambipolar 
diffusion implemented in ZEUS-MP. Therefore, Figure 8 could be a result of case 2) or 3). 
Magnetically subcritical cores are no longer a unique explanation to the B c oc p l J 2 relation. 

In Figures 9a and 9b, we plot the specific angular momentum, j, of each core against 
the radius, R, of the core at t — and 0.65 r ff . The radius R is determined from the 
maximum dimension of the core determined using CLUMPFIND. The majority of cores obey 
the power law relationship, j oc R p , with index p = 3/2, which is plotted for comparison. 
From observations, Goodman et al. (1993) concludes that j scales roughly as R 3 ^ 2 based on 
the velocity gradients of 29 cores in dark clouds. They interpreted this relationship between 
specific angular momentum and core size to imply that cores are in approximate "virial 
equilibrium". With the isothermal equation of state and relatively weak magnetic field 
support in our simulation, the equilibrium in cores is basically between gravitational force 
and rotational angular momentum, producing rotationally supported cores. In their study 
of rotational properties of centrally condensed, turbulent, molecular cloud cores, Burkert & 
Bodenheimer (2000) conclude that the specific angular momentum is related to the radius 
by j oc HR 2 , where f2 is the angular speed of the core. Therefore, j = VLR 2 oc R 3 ^ 2 implies 
f2 oc R^ 1 / 2 . This agrees reasonably well with observations that suggest VL scales roughly 
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R~ 0A (Goodman et al. 1993; Barranco & Goodman 1998). However, we cannot conclude 
just from Figure 9 that the cores are in Keplerian rotation, Q(r) oc r~ a5 , or in angular 
momentum conserved rotation, Q(r) oc r -1 . In addition, there is still a slight scattering of 
data in Figure 9. The rotation could indeed be Keplerian at the inner parts of the core but 
angular momentum conserved in the outer parts. 

In Figure 10a, we plot the relationship between mean density and radius of each core at 
t — (circle) and t = 0.65 rg- (triangle), where mean density < p >= M core /V core , and V core 
is the volume of the core. The result agrees with that of Ballesteros-Paredes & Mac Low 
(2002) (see their Figure 9) that mean density is basically independent of radius. In Figure 
10b, we plot the mass to radius relationship of each core. A power law relationship is seen, 
with M core oc R 2 - 75 from t — 0. Figures 10a and 10b would seem to be inconsistent because 
we would expect M core oc R 3 if cores are spherical. If we have V core oc R 2 - 75 from Figure 10, 
that means either that the geometry of the cores lies between a sphere and a disk, or that 
the cores are a mixture of prolate, triaxial, and oblate objects, as data points are scattered 
around the straight line. We investigate the shapes and orientation of cores in §4.4, and 
confirm this conjecture. The correlations seen in Figures 8 to 10 remain valid for the whole 
simulation even after Truelove et al. (1997) criterion is violated. They are also found at 
lower significance in the 256 3 resolution simulation. 

The number of cores found with CLUMPFIND from t = 0.0 to t — 2.25 basically 
increases with time, as shown in Figure 7. Some cores are destroyed by the supersonic 
turbulence in the simulation. At the same time, some other cores actually merge to form 
more massive cores through gravitation (e.g. Bonnell et al. 1997; Klessen & Burkert 2000; 
Bonnell et al. 2001a,b). Figure 11 shows the merging of two cores, or, to be more specific, 
the accretion of a small clump by a massive clump. The three columns of subplots in Figure 
11 show views along three different axes of the merging sequence. The time interval between 
two rows is 0.05 r s . As the result of core mergers and accretion of surrounding material, the 
final cores possess much larger angular momentum than the initial clumps. This leads to 
the formation of disk-like cores during the simulation. This phenomenon is entirely absent 
in isolated star formation simulations, which consider only a single, gravitationally bound 
core. 



4.3. Core Mass Spectrum Evolution 

It is interesting to examine the effect of core mergers on the core mass spectrum. In 
Figure 12, we plot the evolution of the core mass spectrum from t — to the end of the 
simulation. In each panel, we plot a power-law with slope -1.29 for reference. As discussed 
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in §3.2, at t — 0, the core mass spectrum resembles closely a log normal distribution and 
the slope of the high-mass wing is consistent with the power law with index predicted by 
turbulent fragmentation theory (Padoan & Nordlund 2002). There are some fluctuations of 
the core mass spectra at later time, but the slopes of the high-mass wing remain about the 
same until t ~ Its- The slope becomes shallower at yet later times, but we do not have 
too much confidence in the correctness of the core mass spectra after t = 0.65rfr, as the 
centers of the cores become unresolved. The hydrodynamical simulation by Klessen (2001), 
using SPH and sink particles, also shows shallower slope at the later stages of the simulation, 
which he attributes to the coalescence of cores, but again, fragmentation in the centers of 
cores is suppressed, in this case by the use of sink particles. Even if the trend of increasingly 
shallower slope in core mass spectra is qualitatively correct in our simulation, the rate of 
slope increase may be over-estimated. Our assumed periodic boundary condition means the 
cores have no place to disperse except merging in the later stage. The finite resolution in the 
simulation prevents cores from collapsing as far as real cores would, so that the simulated 
cores have larger cross sections for merging (Paper I). Nevertheless, for a dense cluster, the 
core merging rate may still be high. Recent millimeter wave observations of cloud cores 
also reveal mass spectra of clumps similar to the stellar IMF (e.g Motte, Andre, & Neri 
1998; Testi & Sargent 1998), with slope of the high-mass wing of ~ -1.1 to -1.5, closer to 
Salpeter -1.35 than that derived for gaseous clumps of -0.5 (Blitz 1993), using the form of 
IMF definition in equation (5). Our simulations are in good agreement with these studies. 

Because of the higher resolution of our 512 3 simulation, we have many more cores to 
analyze than in Gammie et al. (2003, see their Figure 4) and we are able to demonstrate 
a core mass distribution similar to the general IMF. Note that not all material within a 
core will finally end up in the protostar. In §5, we will find out that more than half of the 
core mass resides outside the central region. Even material inside the central region may be 
evaporated or ejected away during star formation. In addition, the cores that we observed in 
the simulation may further fragment into binaries or even multiple systems if the resolution 
of the simulation can be increased further. The Truelove et al. (1997) criterion is eventually 
violated at the centers of our cores (see §5) so we cannot follow this evolution. We note 
that, despite the emphasis on artificial fragmentation because of the particular test problem 
chosen by Truelove et al. (1997), regions that are underresolved by their criterion may also 
show too little fragmentation if the geometry is less pathological. Competitive coagulation 
(e.g Silk & Takahashi 1979; Lejeune & Bastien 1986; Murray & Lin 1996; Bonnell et al. 
2001a,b) may reshape the final IMF when gravitation becomes more important. Therefore, 
we cannot reach a definite conclusion that the final stellar IMF is similar to the core mass 
spectrum. 
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4.4. Core shapes and Orientation 

Gammie et al. (2003) recently studied gravitationally collapsed clumps in decaying MHD 
turbulence simulations on 256 3 grids with different magnetic field strengths. In our 512 3 
driven turbulence simulation, va/c s = 1.5, which is close to the value for their run C. Using 
the same definition and symbols of principal axis lengths and ratios as in Gammie et al. 
(2003), we calculated the principal axes of the cores in our simulation at t — 0, 0.65rff , and 
1.31 Tfj and plotted the axis ratios in Figure 13. It is useful to divide the lower triangular 
part of the figure into 3 regions: a prolate region of b/a = to 0.33, a triaxial region of 
b/a = 0.33 to 0.66, and an oblate region of b/a = 0.66 to 1. The result from our higher 
resolution, driven run basically resembles that of Gammie et al. (2003, see their Figure 12). 
The majority of cores are prolate or triaxial in shape, even near the end of our simulation. 
We also see a small number of cores that become oblate late in the simulation, with large axis 
ratios. This implies that some cores become disk-like. Figure 14 shows the angles between 
the shortest core axis and the density weighted mean magnetic field in cores at t — 0.65r ff 
(dot-dash line) and 1.31 r ff (solid line). We see a slight correlation in the angles between 
the shortest core axis and the density weighted local mean magnetic field at t — 0.65r ff and 
an even stronger correlation at 1.31 rg, which is different from the result of Gammie et al. 
(2003). This difference could possibly be a result of higher resolution in our simulation. 



5. Disk-like Cores 

In Figure 13, we see a relatively small number of cores inside the oblate region. The cores 
are too small to study at t — but large enough to study at t — 1.31r ff . Most of these cores 
have a well defined disk-like appearance that could provide some hints about protostellar 
disk formation. Before we begin the discussion of disk-like core, we must emphasize that for 
some cores, the density of the inner cells grows so high that that the local Jeans length is no 
longer resolved in their centers, violating the Truelove et al. (1997) criterion. Unfortunately, 
these cores are also the largest whose physical properties can be studied with reasonable 
statistical significance. Other models using SPH and sink-particles also find disk-like cores 
in the later stages (e.g. Bate, Bonnell, & Bromm 2003). Although our cores are under- 
resolved at t — 1.3lTff, the external disks of cores may still provide some useful information 
on the dynamics of core accretion. With the clear realization that even higher resolution will 
ultimately be needed to properly address the structure of these objects, we proceed with an 
analysis of the data we have in hand, bearing in mind that our 512 3 model is the highest 
resolution driven, self-gravitating MHD simulation currently available. 
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5.1. General Properties of disk-like Cores 

We select 12 cores inside or near the oblate region in Figure 13 with well defined disk- 
like appearance and tabulate their mass and diameters in Table 1. The cores selected are 
shown in the surface density map (Figure 15). Core numbering is arbitrary. The orientation 
of the disk-like cores is uncorrelated to the 3D grid in the simulation despite even though 
they accidentally appear correlated in Figure 15. We have examined the orientation of the 
shortest axis of each disk-like core and they are basically random with respect to the 3D grid. 
The corresponding masses and radii using the high-density scaling mentioned in §2 are also 
shown for easy comparison with observations. The masses of cores range from 0.1 to 2.5 M 
and the diameters of the disks are a few thousand AU. The accretion rates of the cores will 
be discussed later in this section. The 12 cores can be roughly categorized in two groups: 
(1) a cluster of seven cores that remain close to each other in a region with higher mean 
magnetic field, and (2) five cores that are basically isolated and with higher spatial velocity 
(see Table 1). Cores in category 1 have higher mass, central density, magnetic pressure, and 
accretion rates compared to cores in category 2. The mean magnetic field is an order of 
magnitude higher inside the cluster. The size of both grouped and isolated cores appears 
similar. By looking at these two categories of disk-like cores, we can see some manifestation 
of mass segregation in that low-mass cores are preferentially found at large radii from the 
cluster whereas massive cores sink towards the center (e.g. Bonnell & David 1998). 

The mass and sizes of the 12 disk-like cores shown in Table 1, using high-density scaling, 
are in the range of the recent observations of young stellar objects that possess disk-like 
envelopes or circumstellar disks (e.g. Zhang, Hunter, & Sridharan 1998; Fuente et al. 2001; 
Hogerheide 2001; Jayawardhana et al. 2001; Wiseman et al. 2001). Most of these observations 
suggest that the sizes of the circumstellar disks of many young stellar objects are just a few 
hundred AU and the sizes of the disk-like envelopes are a few thousand AU. In addition, 
most of the studied protostars inside the disk-like cores are low or medium mass stars of 
mass < 1 M . The estimated mass of the circumstellar disk is only a small fraction of the 
central star. The more massive objects we find in our simulation are better thought of as 
disk-like envelopes rather than true protostellar disks. 

5.2. Resolution study of disk-like Cores 

In order to better understand how well the cores are being resolved in our simulation, 
we compare them with typical cores formed in other lower resolution simulations. In Figure 
16, we show images of the projected column density of four cores from four simulations at 
the time when the total mass inside cores is ~ 10% of the total box mass. In the 64 3 run 
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( first row of Figure 16), there are four cores found by CLUMPFIND and we cannot identify 
any disk-like core (axis ratio > 2:1). In the 128 3 run, there are 24 cores found, and still no 
disk-like core is identified (second row in Figure 16). In the 256 3 run, there are 38 cores 
determined and four are disk-like (third row in Figure 16). In the 512 3 run, there are 67 
cores and twelve are disk-like. In the last row of Figure 16, we show core 1 from the 512 3 
run. Note that the coordinates shown in Figure 16 are grid zones, not spatial distance. The 
spatial size of the core in 64 3 is the largest but the diameters of cores in 128 3 , 256 3 , and 
512 3 are about the same, indicating the convergence of core size with increasing simulation 
resolution. We can also see that the disk-like appearance of cores becomes more pronounced 
with higher resolution. Figure 17 shows the surface density profiles of the cores in Figure 
16. The curves are offset for ease of comparison. We can see a converging density profile 
of the core with almost uniform density at the center and a power law distribution of the 
disk/envelope with higher resolution, as predicted by other single cloud core simulations. 
Core 1 has the largest spatial diameter and is selected for further study below. 



5.3. Physical Properties of a Disk-like Core 

In Figure 18, we show close-up column density images of core 1 projected along the 
x-, y-, and z-axes. At the end of the simulation, core 1 seems to be an isolated core as 
it is far from other cores. The closest neighbor is about half a box size (> 0.14 pc) away 
from it. In fact, core 1 passes through a dense region at t ~ Its, interacts with several 
cores, and accretes a large amount of material. The disk only appears clearly after this; the 
core remains isolated thereafter. Gravitational interaction among protostars in nascent star 
clusters is argued to tidally truncate or even disrupt accretion disks (e.g. Clarke & Pringle 
1991; Hall, Clarke, & Pringle 1996; Scally & Clarke 2001; Kroupa & Burkert 2001). The 
disk-like cores in this simulation suggest the alternative that disks could survive at least 
some interactions. Core 1 is not massive (~ 0.457 M ) but the extent of the accretion disk 
is the largest (> 5000 AU) of all the cores in the simulation. In the z-direction plot of Figure 
18, we can see some evidence of spiral structure in the disk. 

In Figures 19a-c, the radial profiles of energies, velocities, and surface densities of core 1 
are shown. All the physical quantities in the radial direction are the average values of zones 
binned at the same radius from the center of mass. We can see that the magnetic energy 
remains relatively small compared to the kinetic energy (Figure 19a). The radial velocity 
(dashed line), v ra d, of the disk shown in Figure 19b shows a typical picture of angular 
momentum transfer inside a viscous disk, in which friction tries to spin up the outer part of 
the disk and spin down the inner part (Lynden-Bell & Pringle 1974). However, the largest 
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source of viscosity is likely numerical in our model, so we cannot draw any quantitative 
conclusions from this figure. 

From the choice of sound speed in the simulation, the unit velocity corresponds to 2 km 
s _1 . In Figure 19b, the radial infall rate is < 0.07 in model units, corresponding to < 0.14 
km s^ 1 , on the order of what is observed or modeled (Williams et al. 1999; Myers, Evans, 
& Ohashi 2000). Error bars that plot variances of radial and rotational velocities estimated 
at each radius are also shown. The disk is basically in a rotationally supported state. Note 
that the inward or outward radial velocity of the disk is not the velocity used in estimating 
the accretion rate of material to the core in Table 1, because, as we will show later in this 
section, material mostly accretes onto the core in the direction perpendicular to the disk 
after the disk is well formed. 

The dot-dash line in Figure 19b is the theoretical Keplerian rotation curve (v rot oc r -0 5 ). 
The actual core rotational curve (solid line), v rot , shows typical Keplerian rotation in most 
parts of the disk. The angular momentum conserved rotational curve (dotted line) is also 
shown for comparison. The rotational speed of the core 1 is ~ 0.25-0.35 within the radius 
of 7 zones, which corresponds to ~ 0.13-0.18 rad Kyr -1 in angular velocity. It takes about 
40,000 yr for the core to complete one full rotation. Figure 19c shows the surface density 
profile (solid line) along with the best fitting power law of surface density oc r~ 2A (dash line). 

Observationally, most of the information about true protostellar disks is obtained by 
using the position- velocity (PV) diagram (Richer & Padman 1991). Note that the informa- 
tion on rotation and infall obtained from PV diagrams is subject to large uncertainties if 
the source is poorly resolved. Figure 20 shows the PV diagram of core 1, observed in the 
x-direction along the y-axis (see Figure 19a). The contours in the PV diagram are logarith- 
mic in column density, which is proportional to the observed intensity if the disk is optically 
thin and isothermal. In Figure 20a, the PV diagram is shown at full resolution. The rela- 
tive position (AY) in zones is plotted against the LSR velocity in model units. The shifted 
appearance of velocity in the PV diagram is the result of disk rotation. For qualitative com- 
parison to observations, we convolve the map with a circular Gaussian beam having FWHM 
of three grid zones (Figure 20b). The dashed curve shows the distribution of the rotation ve- 
locity proportional to r _1 , which is characteristic of angular momentum conserving rotation. 
The solid curve shows the distribution of the rotation velocity proportional to r -0 5 , which 
is characteristic of Keplerian rotation. We can see that Keplerian rotation better describes 
most parts of the disk in this core, although the poor resolution of the PV diagram limits 
the strength of this conclusion. If we use this PV diagram to estimate the point mass of the 
central star, we would conclude that the central mass of core 1 is about 0.21 M , using the 
high-density scaling, which is the mass enclosed by the region of radius < 3.4 zones at the 
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center of this core. We find that the mass of the disk outside this central region is more than 
53% of the total core mass (Table 1). An example of an object with properties similar to 
core 1, in high-density scaling, is IRAS 05413-0104 (Wiseman et al. 2001). IRAS 05413-0104 
is a flattened disk-like core of average volume density 10 5 cm -3 , with an envelope of mass ~ 
0.2 M Q , and a disk with diameter of ~ 12,000 AU, roughly twice the size of core 1. 

In Figure 21, we show the velocity and magnetic field components perpendicular and 
parallel to the disk for core 1. The velocity field on the disk plane clearly shows rotation 
(Figure 21b). The edge-on view of the velocity (Figure 21a) shows that material above and 
below the disk is mostly falling onto the core along the magnetic field lines, as expected 
when magnetic field is present. 

The accretion rates of all cores at t — 1.31rg- are listed in Table 1. The accretion rate 
is calculated by integrating the material falling onto the whole area of the disk above and 
below. From Table 1, the accretion rates of the 12 cores are in the range of ~ 10~ 4 -10~ 5 
M yr _1 . If we assume the accretion rate is constant, the ages of the cores are about 15,000 
to 50,000 years old, which is only ~ 0.14-0.46 r ff . Therefore, the accretion rates of cores 
must have dramatically increased near t = 1.31r ff , as expected from a core collapse scenario. 
The rate of accretion onto core 1 is ~ 9.72 xl0~ 6 M Q yr _1 . On average, the cores inside 
the cluster have accretion rate more than twice of those outside. The accretion rates of 
all 12 cores are basically at the high-end of the range of observed envelope infall rates for 
protostellar cores (Ostriker 1998). 

In §4.4, we noted a weak correlation between the core minor axis and the local magnetic 
field direction. Here in Figure 21c, the vertical component of magnetic field lines does show 
a bipolar appearance in core 1, but the magnetic field is bent, due to the movement of the 
core relative to the ambient material. The magnetic field component along the disk also 
shows rotating features, which indicates that the magnetic field is frozen and dragging along 
with the rotating disk (Figure 21d). This magnetic structure is common for all the disk-like 
cores. In Figure 22, we show the 3D structure of the magnetic field in core 6, which is helical 
as a result of core rotation. The angular momenta of cores are expected to be transferred 
to the ambient material by such helical Alfven waves, which in return brake the rotation of 
the cores (e.g. Mouschovias & Paleologou 1979, 1980; Basu & Mouschovias 1994). 

6. Conclusion and Discussion 

Molecular clouds are magnetized and in a state of turbulent motion. It remains an 
open question whether the turbulence is sub-, trans-, or super- Alfvenic (e.g. Mac Low & 
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Klessen 2004). Incompressible hydrodynamic turbulence has been thoroughly studied, but 
compressible MHD turbulence is still a relatively young area, especially as applied to as- 
trophysical problems. In the last decade, numerical studies of compressible turbulence in 
molecular clouds have revealed important results that may directly relate to the early stages 
of star formation. Such simulations still suffer from either lack of magnetic field, lack of self- 
gravity, or insufficient resolution. In this paper, we present the results of our latest studies 
on self-gravitating core formation in a super- Alfvenic turbulent molecular cloud. 

A statistical analysis of the simulation after the turbulence has fully developed shows 
that the power spectra of velocity and magnetic field basically follow Kolmogorov scaling, 
even though the system is compressible and magnetized. This result is consistent with some 
other recent high resolution studies on super-Alfvenic turbulence (e.g. Miiller & Biskamp 
2000; Cho, Lazarian, & Vishniac 2003). The power spectra for density, kinetic, and potential 
energies show shallower power-law indexes. 

We calculated the core mass spectrum in our high-resolution simulation including both 
MHD and self-gravity, and followed its evolution. At least at early times when the model is 
completely resolved, our core mass spectrum is consistent with two main predictions of the 
turbulent fragmentation theory advanced by Padoan & Nordlund (2002): that the slope of 
the velocity spectrum defines the slope of the high-mass wing of the core mass spectrum, 
and that the low mass end flattens and turns over. Padoan & Nordlund (2002) propose that 
this could account for the form of the stellar IMF. At later times in our model, the slope 
of the high-mass wing becomes shallower, possibly because of core coalescence or lack of 
resolution of fragmentation in central regions of cores. Our core encounter rate may well be 
overestimated because of the choice of periodic boundary conditions and grid resolution. 

The simulation that we have presented in this paper is based on initial conditions that 
are both favorable and hostile to the formation of cores. In the case simulated, the molecular 
gas is magnetically supercritical with M/M CT = 8.3. We expect a priori that gravitational 
collapse is unavoidable and will occur eventually. At the same time, the system is subjected 
to a strong turbulent driving force producing flows with rms thermal Mach number M rms = 

10. The turbulence is strong enough to formally support the entire region, so it could prevent 
or delay the collapse of cores or even destroy some cores in the process of collapsing. In Paper 

11, it was shown that in a simulation on a 256 3 grid, with the same initial conditions, the 
magnetic and turbulent support could not prevent local collapse of cores in the molecular 
cloud. Cores form almost immediately after the gravitational force is turned on. One of 
the motivations for our higher resolution simulation — twice the spatial resolution and eight 
times the mass resolution of the 256 3 model of Paper II — is to determine whether stronger 
Alfven wave support may occur when shorter wavelengths are resolved in a simulation using 
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the same initial conditions. We find that gravitational collapse still wins over magnetic and 
turbulent support. The 512 3 simulation lasts about 2/3 of a full free-fall time, after the 
gravitation is turned on, before the first cell in the computation violates the Jeans criterion 
Truelove et al. (1997) that the local Jeans length be resolved, a factor of four longer than 
the 256 3 simulation. A resolution study shows convergent trends with increasing simulation 
resolution of both the fraction of the box mass in bound cores (Fig. 5) and the core size 
(Fig. 17). 

Our simulation began with the entire box magnetically supercritical. However, it is 
often argued that gravitational fragmentation will produce smaller and smaller regions that 
eventually become subcritical (Mestel & Spitzer 1956). We found that gravitationally bound 
cores formed and collapsed until they were rotationally supported (Fig. 9), while still having 
masses at least an order of magnitude over the critical mass for magnetic support (Figs. 
8a-b). We further found that the central magnetic field strength depends on the central 

1/2 

density as B c oc p c , suggesting that observations of such relations do not necessarily point 
to subcritical cores. 

Most of the cores appear to be prolate or triaxial in shape, consistent with simulation 
results from Gammie et al. (2003). However the correlation of the shortest axis to the 
density-weighted local magnetic field of the cores is found to be stronger in our simulation. 

In our simulation, we observe interaction between cores during close encounters. The 
continuing accretion of large amounts of material, and merging of cores results in cores with 
high specific angular momentum that helps the formation of disk-like structure. The final 
outcome is highly unpredictable as both magnetic and turbulent forces are present during 
the gravitational collapse process. In our simulation, we have identified twelve disk-like cores 
with axis ratio > 2:1 and large enough for statistical analysis. The sizes and mass of the 
twelve cores are in the range of young protostellar objects observed with flattened disk-like 
structure (Table 1). All the disk-like cores are rotationally supported. Magnetic pressure 
support is relatively unimportant for these cores, perhaps due to the large M/M cr ratio that 
we chose as an initial condition. The surface density of most of the disk-like cores has a 
power law distribution, as observed in real cores, and as predicted in isolated core collapse 
simulations. Instead of being tidally truncated or disrupted, the core disks survive and 
flourish among strong interactions. 

Most of the disk-like cores appear to be in Keplerian rotation from their PV diagram. 
Observationally, many protostellar disks are found to be in Keplerian rotation (e.g. Wiseman 
et al. 2001, and references therein). There are still some observed protostellar disks found 
to be in non-Keplerian rotation (e.g. Myers, Evans, & Ohashi 2000). Wiseman et al. (2001) 
points out that the non-Keplerian rotation curve reflects a mass distribution that is not yet 
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dominated by a central condensation. Also, a system could be in Keplerian rotation but a 
high optical depth would result in velocity maps that are dominated by the motion in the 
outer layers along the line of sight toward the center of the core. For the twelve rotationally 
supported disk-like cores, a large amount of matter is still accreting perpendicular to the 
disks (Figure 21a); converted to astrophysical units this corresponds to substantial accretion 
rates (see Table 1). When the mass of the disk keeps increasing, the gravitational stability 
of the rotational disk may be affected. Unfortunately, because of the limited resolution, we 
cannot follow the core collapse long enough to see consequences beyond the first hint of spiral 
structure in the disk. 

There is as yet no direct detection of disks around high-mass (> 1 M Q ) protostars. 
Fuente et al. (2001) suggests that for stellar mass > 5 M , destruction mechanisms like 
radiation pressure on dust or ionization by stellar ultraviolet emission may be responsible 
for the rapid erosion of the outer disk. However, these mechanisms would only operate at 
the later stages of protostellar evolution. Zhang, Hunter, & Sridharan (1998) point out 
that the lack of direct detections of disks around high-mass protostars may simply be due 
to the greater distances to high mass star formation regions, as well as their clustering, 
which further complicates the kinematics. Therefore, the existence of disks around high- 
mass protostars at the very early stage of core collapse remains uncertain. It is important to 
note that in a number of objects studied at high resolution with radio and millimeter-wave 
interferometry, flattened structures and Keplerian motions argue strongly for the existence of 
rotationally supported disks (e.g. Sargent & Beckwith 1991; Koerner, Sargent, & Beckwith 
1993; Dutrey 1996). Improvements in the resolution of both observations and simulations 
on the early stage of core collapse in molecular clouds are vital to understand the details of 
star formation inside magnetized turbulent molecular clouds. 
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Fig. I. — (a) Power spectra of density p (solid line), velocity v (dashed line), potential 
energy Up (dotted line), kinetic energy Uk (dash-dot-dot-dotted line), and magnetic energy 
Up (dash-dotted line) of the 512 3 simulation. The thin straight line is the Kolmogorov power 
spectrum, which fits well to the inertial range of velocity and magnetic energy power spectra. 
The long dashed line is the driving spectrum, (b) Velocity power spectra of 256 3 (circle) and 
512 3 (triangle) simulations of the turbulence, compensated by /c 5//3 , at t — when gravity is 
turned on. The horizontal straight line is the Kolmogorov power spectrum. 
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Fig. 2. — Density PDF and log-normal fit in 256 3 (open square and dashed line) and 512 3 
(open circle and solide line) simulations at t — 0. 




Fig. 3. — Column density of the 512 3 simulation projected along the x-, y-, and z-directions. 
All images are plotted using the same gray scale. A log column density of -0.602, corresponds 
to gas at the mean density times the box length, (a - c) At t — 0, when gravitation is switched 
on, filamentary structure dominates the appearance, (d - f) At t — 0.65rff, filaments break 
into pieces and higher density clumps slowly increase their mass by accreting material around 
them, (g - i) At t — 1.31r ff , many cores with high column-density contrast are formed. 
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Fig. 4. — 3D volume rendering of the logarithm of density (left) and the magnetic pressure 
(right) of the turbulence box at t — ISItq. The whole simulation is also available as an 
mpeg animation in system time unit. 
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Fig. 5. — The mass spectrum of cores (dot-dashed line) and clumps (dashed line), in the 
256 3 simulation and the mass spectrum of cores (solid line) in the 512 3 simulation at t — 0. 
The spectra are normalized by the Jeans mass, Mj(p). The straight lines are the power law 
with index -1.29 predicted by Padoan & Nordlund (2002) for our velocity power spectrum. 
High-mass wings of all 3 spectra are consistent with the predicted power law. 
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Fig. 6. — Comparison of the mass in gravitationally bound cores for runs driven at k — 
1-2 with varying resolution 64 3 (triangle), 128 3 (square), 256 3 (plus), and 512 3 (circle). 
M* denotes the sum of masses found in all cores, in units of box mass, determined by the 
modified CLUMPFIND (Williams, De Geus, & Blitz 1994, ; see also Paper I). Collapse rates 
vary but collapse occurs in all cases. 
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Fig. 7. — Core number counted from the simulation using the modified CLUMPFIND. The 
decrease in core number is caused by the destruction of cores by supersonic turbulence or 
merging of cores. 
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Fig. 8. — Evolution of magnetic field and density of cores. Ratio of core mass to critical 
mass for magnetostatic support M/M cr is shown for cores at (a) t — 0, and (b) t = 0.65r ff 
in x- (circle), y- (triangle), and z-direction (square). Cores are magnetically supercritical. 
Central magnetic field strength B c vs. central density of cores p c is shown at (c) t — 0, and 
(d) t = 0.65r ff . The straight line is a power law of index 0.5. 
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Fig. 9. — Evolution of specific angular momentum j vs. radius R of cores at (a) t — 0, and 
(b) t = 0.65rff. The straight line is a power law of index 3/2. The data distribution matches 
well with observations of dark cloud cores (e.g. Goodman et al. 1993). 
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Fig. 10. — (a) Average density < p > vs. core radius R at t — (open circles) and t = 0.65 
Tgf (open triangle). Basically, core radius is unrelated to core average density (Ballesteros- 
Paredes & Mac Low 2002). (b) Core radius R vs. core mass M core at t — (open circles) 
and t = 0.65r ff (open triangle). The straight line is a power law of index 2.75. All quantities 
are in scaled system units. 
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Fig. 11. — Different views of a core accretion event along the x-, y-, and z-directions. The 
time interval between two rows of the figures is 0.05 Tg. Core merging and accretion of large 
material clump are an important process in forming disk-like cores. 
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Fig. 12. — Evolution of the core mass spectrum, normalized by the Jeans mass Mj(p ), 
in the 512 3 simulation. The straight line is the power law with index -1.29 predicted by 
Padoan & Nordlund (2002) for (5 — -5/3 (see text). The core mass spectra generally show 
a "universal" IMF appearance (Kroupa 2002) in the simulation with a clear turn over at 
the low mass region. The slope of the high-mass wing of the spectra matches well with the 
turbulent fragmentation prediction and remain about the same until t = 0.49rff. Later the 
slope becomes shallower as the result of core coalescence. 
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Fig. 13. — The distribution of core axis ratios at t — (square), 0.65 (circle), and 1.31 r ff 
(triangle). Cores near the diagonal are prolate and the cores near b/a — 1 are oblate. The 
region is divided into 3 parts and the cores inside the middle region are classified as triaxial 
(Gammie et al. 2003). The majority of cores are prolate or triaxial in shape. The principal 
axes of cores are defined such that a > b > c. 
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Fig. 14. — The correlation of angles between the shortest body axis (e c )and the density 
weighted mean magnetic field (Bi) in the cores at t = 0.65rg (dot-dashed line) and 1.31 rg 
(solid line). We observe stronger correlation of the angles in this simulation than that in the 
simulation by (Gammie et al. 2003) with similar initial weak field. 




Fig. 15. — Column density map highlighting the location of 12 disk-like cores at t 
The properties of core 1 is analyzed in detail in §5.3. 



= 1.3lTff. 
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Fig. 16. — The effect of numerical resolution on core structure when total 10% of material 
is trapped inside cores: 64 3 (1st row), 128 3 (2nd row), 256 3 (3rd row), and 512 3 (4th row). 
Grey scale images of column density, projected along x-, y-, and z-axes and the axes denote 
cell number. Disk-like core structure barely seen in 256 3 resolution and becomes clear in the 
512 3 simulation. 
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Fig. 17. — Surface density of the cores in four simulations with resolutions 64 3 , 128 3 , 256 3 , 
and 512 3 , as shown in Figure 16. The surface densities of the cores in 64 3 , 128 3 , and 256 3 
resolution simulations are scaled by 10 4 , 10 3 , and 10, respectively for comparison. With 
higher resolution, the core surface density profile converges to a uniform density central 
region and a power law disk. 
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Fig. 18. — Logarithmic column density images of Core 1 in 3 directions in gray scale and 
contours. The contour interval is 0.265. A spiral structure of the accretion disk is barely 
visible in the z-direction. 
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Fig. 19. — Radial profiles of spherically averaged properties of core 1 at t — 1.31rg. (a) 
Kinetic (solid line), magnetic (dashed line), and potential (dot-dashed line) energy distribu- 
tion along the radius, (b) Rotational (solid line) and radial (dashed line) velocity profiles. 
The theoretical rotational curves for a Keplerian rotation (dot-dashed line) and angular mo- 
mentum conservation rotation (dotted line) are also shown. Keplerian rotation seems to fit 
better to the data, (d) Surface density of the core shows a power law outer profile. The 
slope of the power law fit (dot-dashed line) is -3. 
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Fig. 20. — Position- velocity diagrams of the Core 1. (a) PV diagram before convolution. 
The contours are the logarithm of column density, (b) PV diagram after convolving (a) with 
a circular Gaussian beam with a F WHM of 3 zones. The solid curve indicates the Keplerian 
rotation and dotted curve indicates the angular momentum conservation rotation. 
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Fig. 21. — Vector plots of velocity and magnetic field in Core 1. (a) Side view along the x- 
direction shows material accretion along the polar direction, (b) Plan view along z-direction 
shows disk rotation, (c) Side view along the x-direction shows a distorted bipolar structure 
of the magnetic field, (d) Plan view along z-direction shows dragging of magnetic field by 
the core rotation. Surface density images are shown in gray scale. 
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Fig. 22. — The 3D view of the helical magnetic field through core 6 caused by core rotation. 
The solid rendered objects are density isosurfaces. 
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Table 1. Physical properties of 12 disk-like cores. 



Core 


Category 


Mass a 


Diameter a 


Accretion Rate a 


Mass b 


Diameter* 3 


Accretion Rate b 












(M ) 


(AU) 


(M Q yr- 1 ) 


1 


1 


0.0035 


0.0977 


5.25X10 2 


0.457 


5642 


9.72xl0~ 6 


2 


2 


0.0114 


0.0586 


2.18X10 3 


1.489 


3384 


4.04xl0~ 5 


3 


1 


0.0028 


0.0547 


3.94X10 2 


0.366 


3159 


7.30xl0~ 6 


4 


1 


0.0033 


0.0508 


2.07X10 3 


0.431 


2934 


3.84xl0~ 5 


5 


2 


0.0192 


0.0625 


3.91 xlO 3 


2.508 


3610 


7.25 xl0~ 5 


6 


2 


0.0121 


0.0547 


2.73x10 s 


1.580 


3159 


5.06xl0~ 5 


7 


2 


0.0015 


0.0430 


1.26X10 3 


0.196 


2483 


2.34xl0~ 5 


8 


1 


0.0017 


0.0586 


1.53X10 3 


0.222 


3384 


2.83X10- 5 


9 


1 


0.0011 


0.0313 


1.05X10 3 


0.144 


1808 


1.95xl0~ 5 


10 


2 


0.0064 


0.0586 


1.90X10 3 


0.836 


3384 


3.53xl0~ 5 


11 


2 


0.0048 


0.0781 


2.22X10 3 


0.627 


4511 


4.12xl0~ 5 


12 


2 


0.0138 


0.0195 


5.59X10 3 


1.802 


1126 


1.04xl0~ 4 



a Units in simulation scaling 
b Units in high density scaling 

Note. — Cores in category 1 belong to a group with density, magnetic pressure, and accretion rate 
generally higher than those of category 2 cores found in open space. 



